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Introduction 

Dendritic growth is important because it is the basic microstructural pattern in solidified metals. The 
pattern selected during solidification of a pure material depends on the existing thermal field during 
freezing. Once this pattern is set, it is difficult to change in the solid state without substantial effort, e.g., 
through mechanical deformation and heat treatment. 

The evolution of dendritic microstructures is reasonably well understood for pure materials growing 
into undercooled melts under purely diffusive conditions, i.e., when fluid flow is absent. The tip of the 
dendrite approximates a paraboloid of revolution. Ivantsov [1] detennined the solution for the thermal 
field surrounding an isothermal dendrite, neglecting surface tension, and determined a relation between 
undercooling and the (constant) tip velocity and tip. The function is a combination of exponentials and 
error integrals, whose form is well known and not important for the current discussion. Since the tip radius 
and velocity appear only as a product, the diffusion solution does not uniquely determine the shape. 

The pattern selection problem is resolved by including surface tension and its anisotropy in the boundary 
condition for the temperature of the dendrite, and relaxing the assumption that the shape is known. [2,3] 
This body of theory is known as “microscopic solvability.” In 2-D, it gives a correction to the shape of the 
interface near the tip, converging to the Ivantsov solution far away from the tip where curvature becomes 
negligible. The theory provides a second relation for the dendrite tip velocity and radius. The problem is 
somewhat more complicated in 3-D, where corrections to the Ivantsov shape are large. See Brener [4] for 
a treatment of this problem. 

Recent numerical calculations using the phase field method, in two dimensions and at high undercooling, 
agree very well with the predictions of microscopic solvability theory. [5,6] At low undercooling, the 
computed selection constant still agrees with the theory, but the tip radius and tip velocity differ. Provatas 
et al [6] showed that this discrepancy arises because the transport solution at low undercooling does not 
satisfy the premise, assumed in solvability theory, that the dendrite is a single, isolated branch, independent 
of its neighbors. This may explain the difference between the predictions of microscopic solvability 
theory and experimental observations, which are invariably made at low undercooling. 

The microstructure is significantly altered by the presence of flow during solidification. [7,8] This problem 
is important because flow induced by buoyancy, residual pouring currents or forced flow is nearly always 
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present in castings unless great pains are taken to avoid it, for example by performing experiments in the 
reduced gravity of outer space. Saville and Beaghton [9] extended the theory of Ivantsov to consider the 
case of an isolated paraboloid of revolution, with zero surface tension, growing in a shape preserving 
way at constant velocity, with a uniform forced flow from infinity aligned parallel to the growth axis. 
The characteristic parameters for this problem are such that the flow is dominated by viscous forces, but 
the transport is dominated by advection. The flow is therefore well approximated by Oseen’s equation, 
and Saville and Beaghton presented a solution where the tip velocity and radius are still not uniquely 
determined by the transport solution, and a selection criterion is still required. 

Bouissou and Pelce [10] examined the stability of a 2-D dendrite under the assumptions of Saville and 
Beaghton’s analysis, and determined the scaling of the selection constant with flow. The applicability of 
this theory is not yet confirmed by either computations or experiments. A recent study by Beckermann 
et al [11] found a weak dependence of the selection constant on flow, but neither confirmed nor negated 
the theory. The existing experimental evidence is also inconclusive, and perhaps contradictory. Bouissou 
et al [12] examined the flow of a dilute pivalic acid-ethanol (PVA-EtOH) solution past a solidifying PVA 
dendrite, and found a trend similar to that predicted by their theory. On the other hand, in a series of 
experiments using pure succinonitrile (SCN), Lee et al [ 1 3] found that the selection constant increased with 
flow velocity. These discrepancies are currently unexplained, and in this paper we perform simulations 
intended to shed some light on this issue. We use the phase field method, which we now describe briefly 
before proceeding to a more detailed description of our simulations. 

The phase field method has become the method of choice for simulating dendritic growth, because of its 
ability to handle the complex evolving shape of the dendrite. [14,15,16,17,18] Because this method has 
been described in detail in numerous articles, we present here just a brief sketch. The dendrite growth 
problem follows the evolution of the solid from the liquid. In the so-called sharp interface problem, the 
liquid-solid interface must satisfy two boundary conditions: thermodynamic equilibrium, which requires 
that the interface temperature satisfy the Gibbs-Thomson condition including crystalline anisotropy 

Solution of the sharp interface problem for dendritic growth is difficult, because it requires that the moving 
boundary be tracked explicitly. In the phase field method, a continuous order parameter is introduced, 
such that -1 corresponds to the liquid, +1 corresponds to the solid, and the set of locations where the 
parameter is 0 corresponds to the interface between the two phases. The liquid-solid interface is diffuse, 
with nominal thickness W Q . For a more detailed discussion of the method, the specific models used, and 
a description of convergence of the phase field model to the sharp interface, see Karma and Rappel [5], 
Beckermann et al [11] and Jeong et al [19]. 

The phase field method introduces an artificial finite width for the interface. Karma and Rappel [5] 
demonstrated convergence of the phase field model to the sharp interface model in the limit where the 
interface width is much less than the diffusion length, and gave specific formulas for relating the parameters 
in the two problems. Grid convergence of the numerical solution requires that the grid space\ing be much 
less than the interface width and also that the domain size be much larger than the diffusion length. See 
Karma and Rappel [5] for further details. The important result is that these two requirements force us to 
choose the domain to be 1000 or more times the grid spacing. Thus, a uniform mesh would require on the 
order of 10 6 grid points in 2D, and 10 9 in 3D. 
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We finesse the length scale problem by solving the equations on an adaptive grid, providing high resolution 
near the interface, and a coarser mesh further away to resolve the diffusion field. The methods have been 
described elsewhere [6,19], and so we provide only a limited discussion in this article, in the next section. 
Plapp and Kanna [20] developed a different approach, where they solve the phase field model on a regular 
inner grid, and resolve the outer diffusion field using a Monte-Carlo method, with the two solutions 
matched at an artificial boundary. Both methods are effective for the diffusion problem, but the extension 
of the latter method to cases where fluid flow is present is not obvious. Accordingly, we focus the rest of 
our discussion on the solution of the solidification problem with melt convection using our adaptive grid 
techniques. 

It is important to understand that in order to compare our numerical calculations with experiments, we must 
perform simulations in three dimensions, at low undercooling, and for materials of low anisotropy. Each 
one of these choices makes the problem computationally more difficult, and the combination represents a 
formidable challenge. We have developed a code which solves the 3D Navier-Stokes, energy and phase 
field equations on an adaptive grid in parallel. [19] This is essential for us to take up the problems we 
have described. 

In the next section, we give an abbreviated description of the phase-field model, and its numerical 
implementation. We then describe the results of simulations relevant to the experiments of Lee et al [13] 
and Bouissou et al [12]. 

Methods 

The phase field method has been extended to include fluid flow in the melt along with solidification. 
[21,22,11,19] The essential new features, in addition to solving the Navier-Stokes equations in the fluid 
phase, are the formulation of the model such that the velocity is extinguished in the solid phase, and 
ensuring that the interfacial shear stress is correctly represented. Tonhardt et al [21] and Juric [22] use an 
enhanced viscosity approach to extinguish the velocity in the solid phase. In this approach, the viscosity 
changes rapidly from its value in the melt to a very large value representing the solid, as the phase field 
goes from - 1 to 1 . This method is effective, but convergence can be very sensitive to the relative magnitude 
of the viscosity change and the interface thickness. 

We adopt the method of Beckermann et al [11], who introduced a mixture formulation for the continuity, 
Navier-Stokes and energy equations. An interfacial stress term is added to the Navier-Stokes to ensure 
that the shear stress is correct at the interface. See Beckermann et al [1 1] for further details. 

The important physical parameters in the simulations are the thermal diffusivity, capillary length and surface 
tension anisotropy. The phase field model introduces the additional non-physical parameters of the grid 
spacing, interface width and domain size. Karma and Rappel [5] presented an asymptotic analysis to show the 
relations between the parameters such that the phase field model converges to the sharp interface model. 

We solve the 3D flow equations using the Semi-Implicit Approximate Projection Method (SIAPM)[23], a 
predictor-corrector method which can solve these equations effectively, especially for large 3D problems, 
using relatively small amounts of memory. For a detailed discussion of the algorithm, the reader is 
referred to the original paper [23], and for a detailed description of the parallel adaptive finite element 
implementation, see Jeong et al [19]. 
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The equations are solved in a segregated fashion, using an implicit time-stepping scheme, except for the 
phase-field equation, which we solve explicitly. The grid is adapted so that it is finest in the vicinity of the 
interface, where the grid spacing must resolve the interface width W , to a much coarser mesh far away, 
where all of the primitive fields vary slowly. The typical ratio of largest to smallest element size in our 
calculations is 512, but there is no inherent limit. Element refinement and fusion is done based on the 
value of an error estimator, computed from the solution within each element. The details of the procedure 
are given in Jeong et al [19], and are not repeated here. The grid was adapted whenever the dendrite tip 
moved a predetermined distance, typically 0.41^. 

The basic problem which we solve places a small spherical seed at the centroid of a cube, with the axes 
of crystalline symmetry aligned with the cube edges. Earlier work has demonstrated that when the grid 
spacings are chosen as described earlier, grid anisotropy is minimal, and the seed orientation is merely a 
matter of convenience. Flow is introduced on the front face of the cube as a uniform velocity. The lateral 
faces are all assumed to be planes of symmetry, and the exit face has zero shear stress. The temperature 
on all faces of the computational domain is set to the undercooling, and we always use a domain of edge 
length sufficiently large that the thermal and velocity boundary conditions do not affect the result. 

A convenient measure of the results is to track the evolution of the tip velocity and tip radius in the 
various directions. Simulations in both 2D [11] and 3D [19] showed that the tip velocity is increased on 
the leading arm, decreased on the trailing arm, and essentially unchanged on the transverse arms. We 
introduced a few “tricks” to speed up the calculations, and these are described now, using the leading and 
trailing tip velocities as indicators of the validity of these techniques. 

We use a different time step size to solve the various equations. The phase-field and temperature equations 
are solved using the minimum time step and the fluid flow equations are solved using a larger time step. 
The physical basis for this choice is that the interface moves very slowly relative to the fluid, and thus the 
motion of the interface does not strongly affect the flow field. The results were indistinguishable for ratios 
up to about 10 for relatively high speed growth, and up to a ratio of 25 for very low speed growth such as 
in SCN at low undercooling. This provides a substantial computational savings compared to using a single 
time step for all of the fields. In order to further reduce computation time, we also use a technique where 
we start the computations with flow using a seed that had been grown to small, but finite size without flow. 
The tip velocity for a dendrite started from the precursor seed quickly recovers to the solution where flow 
was present from the beginning, and the results are indistinguishable. The dendrite tip radius shows the 
same trend. 

Results 

The simulations we report in this section were run to examine the role of fluid flow on dendrite growth at 
low undercooling, for comparison with experiments. We present results for the cases shown in Table 1. 
The value of anisotropy equal to 0.0055 corresponds to SCN, while anisotropy equal to 0.025 corresponds 
to PVA. We note that these values have not been corrected for grid anisotropy. The results presented below 
should nevertheless be viewed as approximate in this sense. 
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Table 1 : Parameters for simulations reported in this section 


Case 

Anisotropy 

Undercooling 

U (cm/s) 

la 

0.0055 

0.04 

0 

b 



0.5 

c 



1.0 

2a 

0.025 

0.05 

0 

b 



0.12 

c 



0.24 

d 



0.60 


The simulations were run on a variety of machines: Sun Ultra-2/200 MHz, IBM RS/6000 7043 Model 
260, and SGI Origin 2000 multiprocessor machines. Comparisons of cpu time for the various cases are 
therefore not meaningful. The parallel runs on the SGI Origin 2000 were usually made as a series of 
restart/batch jobs in order to optimize throughput in the queing system. Statistics for one run (Case lc), 
which should be considered representative, are that the run took 2,496 cpu hours on 32 processors, 101 
clock hours (average parallel efficiency of 77%), and used about 2 Gb of memory distributed over the 32 
processors. 

Succinonitrile 

With flow, the tip velocity increases dramatically, and tip radius decreases dramatically, immediately as 
the flow is started. The simulations do not appear to have reached a final steady state, but the temporal 
variation is small at the end of the simulations. (The simulations were terminated when the mesh size 
grew to about 450,000 nodes). The selection constant reaches its final value much more quickly than either 
the tip velocity or radius, which is also characteristic of calculations in the absence of flow. We find that 
the selection constant decreases slightly with flow. 

Table 2 compares the computed values with those obtained experimentally by Lee et al. [13] The results 
are presented as ratios of the final value in the simulations to the corresponding value without flow. Note 
that our simulations are not fully converged. Whereas the ratios of tip velocities and tip radii are quite 
similar in both the simulations and the experiments, the trends in the selection constant are opposite, i.e. 
we find it to be a weakly decreasing function of the far-field velocity, while Lee, et al. found it to be an 
increasing function of flow velocity. 

Table 2: Comparison of computed and experimental tip parameters for succinonitrile. 



Flow 

Tip velocity 

Tip radius 

Selection 


velocity 

(cm/s) 

ratio 

ratio 

constant 

ratio 

Expt 

1.0 

2.2 

0.6 

1.4 

Calc 


2.6 

0.72 

0.8 

Expt 

0.5 

1.7 

0.71 

1.2 

Calc 


2.1 

0.71 

0.95 
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Finally, it is in instructive to compare the results of the calculations with the solution of Beaghton and 
Saville. Although there is no unique state determined by their theory, the computed solution might be 
expected to follow their transport solution. When the tip radius is computed using the exact computed 
shape of the dendrite, the results are quite far from the Oseen-Ivantsov solution. Beckermann et al. 
[11] showed that the tip radii should be corrected to compute a smooth paraboloidal shape near the tip 
for comparison to the transport solution. For this case, however, because the tip is so broad, the shape 
correction has almost no effect. 

While the trends are similar, none of the curves matches very well. This is important, because it shows 
that neither the experiments nor the calculations match the transport solution for an isolated paraboloid 
of revolution. We believe that the explanation for this is that neither the computed nor experimental 
dendrites are well approximated as paraboloids, and therefore the temperature field should not be expected 
to match the transport solution. It is also possible that the transverse arms affect the flow, an effect not 
accounted for at all in the Oseen-Ivantsov theory. We observed this phenomenon in diffusive growth at 
low undercooling [24]. 

PVA 

We also considered several cases at higher anisotropy of 0.025, corresponding to PVA. At high 
undercooling (0.55), our results are similar to those reported by Beckermann et al, whose simulations 
were done in 2D, but were otherwise similar to our 3D analyses. The results compare reasonably well to 
the Oseen-Ivantsov. In this case, where the dendrite tip is much sharper than it was for SCN, the parabolic 
tip curvature correction brings the results into agreement with the Oseen-Ivantsov transport solution. This 
occurs because at high values of anisotropy and undercooling the computed shape tends to look much more 
like the isolated paraboloid in the analytical solution, than it does for smaller values of the parameters. 

We note that we cannot compare these calculations directly with the experiments of Bouissou et al, for 
several reasons. Their experiments were performed in an alloy system, and at present we compute only for 
pure materials. Further, their experiments were performed by growing dendrites in a narrow gap, much 
thinner than the boundary layers in that direction for a freely growing dendrite. We will examine this case 
in a later paper. 

For PVA at low undercooling, similar trends are observed for PVA as in SCN, i.e. the dendrite tip velocity 
increases and tip radius decreases as flow velocity increases. Once again, the selection constant is a 
weakly decreasing function of flow velocity, which reaches a steady value much sooner than either tip 
radius or velocity. The isotherms are once again advected by the flow, but the tip shape is much smoother 
for PVA than for SCN. There are still systematic deviations from the Oseen-Ivanstov solution. 

Discussion 

The most important observation from our computations is that significant discrepancies exist between 
the analytical theory of Saville and Beaghton for an isolated branchless dendrite, our calculations and 
the available experimental results. We believe that the source of these differences lies in the difference 
between the real shape of a growing dendrite, and the assumed branchless, isolated paraboloid. 

Much has been made about the trend in the selection constant with increasing imposed flow. We find a 
ver/ weak negative dependence of on flow velocity, i.e., the selection constant decreases slightly as flow 


velocity increases. We note that even though the selection constant is nearly constant, the tip velocity 
increases substantially, and the tip radius decreases significantly. Within the precision of our calculations, 
the selection constant is essentially constant. This result is consistent with the theory of Bouissou and 
Pelce. 

There were significant deviations between the theory of Saville and Beaghton and both our calculations 
and the experiments of Lee et al [13]. We believe that the difference arises from the idealization of the 
dendrite as an isolated paraboloidal shape. For low Reynolds numbers, we expect to find a viscous 
boundary layer ahead of any obstacle, including sidebranches and transverse arms. For the example cited 
in the preceding paragraph, the boundary layer is larger than the distance to the transverse arms in our 
computations, and probably also in the experiments. Thus, if the flow field is significantly different from 
the Oseen-Ivantsov solution, one cannot expect the results to match the theory. The flow reacts to the 
entire body of the dendrite, not just the tip. We note that this effect would be even more pronounced in 
two dimensions, where flow around the dendrite tip cannot occur. 

Our conclusion from this is that both the experiments and the calculations are very difficult to perform in 
such a way as to satisfy the idealized conditions assumed in the theory, and this clouds our ability to make 
precise statements about the theory of tip selection under conditions of forced flow. 

Conclusion 

We have developed a three-dimensional, adaptive, parallel finite element code to examine solidification 
of pure materials under conditions of forced flow. We have examined the effect of undercooling, surface 
tension anisotropy and imposed flow velocity on the growth. The flow significantly alters the growth 
process, producing dendrites that grow faster, and with greater tip curvature, into the flow. The selection 
constant decreases slightly with flow velocity in our calculations. 

The results of the calculations agree well with the transport solution of Saville and Beaghton [9] at high 
undercooling and high anisotropy. At low undercooling, significant deviations are found. We attribute 
this difference to the influence of other parts of the dendrite, removed from the tip, on the flow field. 
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